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1. Introduction 


Chebyshev Spectral Tau methods are a useful technique for solving ordinary dif- 
ferential equations. Numerical programs using this method axe often considerably 
faster with greater accuracy than other standard techniques such as finite differ- 
encing. In my attempt to understand Chebyshev polynomials and their applica- 
tions, I was unable to find one definitive book or article which took me from only 
a general knowledge of calculus and linear algebra to being able to confidently 
apply these new techniques. This tutorial is an attempt to fill that need. Several 
examples are worked out to demonstrate the usefulness and occasional difficulty 
of their application. 

This series has been broken down into four major parts. Section one gives 
a brief explanation and example of when and why Chebyshev Polynomials are 
so useful. The second section is a detailed account of the definitions, properties 
and algebra that surrounds Chebyshev Polynomials. Because section two is so 
detailed, the reader may read only the first few and last few paragraphs to get 
an overview of its material. Section three is also a detailed piece. It deals with 
the recurrence relationship between functions and Chebyshev Polynomials. These 
recurrence relationships are very important in the implementation of Chebyshev 
Polynomials to spectral methods. None the less, most of the common recurrence 
relationships are tabulated in the appendix of Gottlieb and Orszag’s book [4], For 
this reason, the hurried reader may just browse over the third section. Section 
four gives several examples of how to implement Chebyshev Polynomials with 
spectral tau methods. All the examples are linear differential equations and most 
involve finding the eigenvalues of these equations. 

It can be shown that the coefficient on the x n term of the Chebyshev polyno- 
mial is 2 n_1 . This is left as an exercise in Exercise 2.3. In this tutorial we will 
take some function, truncate it and expand it in terms of a Chebyshev Polynomial. 
When a function is expanded in terms of some other polynomial, the coefficients 
of the expansion are roughly inversely proportional to the coefficients of the poly- 
nomial. Therefore, a function expanded in terms of a Chebyshev Polynomial will 
have coefficients roughly proportional to 2 1_n . This gives a geometric rate of con- 
vergence as the number of terms in the expansion, n, is increased. The following 
example 1 is taken from an excellent source for spectral methods , Spectral Meth- 
ods in Fluid Dynamics [2]. The example deals with a linear heat equation with 

1 Reproduced with authors permission 
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homogeneous Dirichlet boundary conditions, over the interval of [—1,1], 


du d 2 u 

(1.1) 

dt dx 2 

u(l,t) = 0 

(1.2) 

u(— l,t) = 0 

(1.3) 

u(x,0) = sinra 

(1.4) 


A Chebyshev collocation method and a second order difference equation axe used 
to solve for u. The convergence of the two compared to the number of terms is 
given in the table below. 


Table 1.1 Maximum error for the one-dimensional heat equation 


N 

Chebyshev collocation 

Second Order finite difference 

8 

4.58 * 10- 4 

6.44 * 10" 1 

10 

8.25 * 10' 6 

3.59 * 10- 1 

12 

1.01 * 10" 7 

2.50 * lO -1 

14 

l.io* nr 9 

1.74* 10" 1 

16 

2.09 * 10~ u 

1.35* 10- 1 


Looking at the errors between the two methods, it becomes rather obvious that 
the Chebyshev method converges significantly faster. This faster convergence rate 
would produce enormous savings in computational time. 

2. Chebyshev Polynomials 

2.1. Definition 

The Chebyshev polynomial is defined as 


T n (x) 

= cos [n arccos(x)] 

(2.1) 

-1 

< X < 1 

(2.2) 

X 

= cos 9 

(2.3) 


2.2. Minima and Zeroes 

The minima are the values in which the Chebyshev polynomial is equal to zero. 


T n (x ) = 0 = cos(n9) 
cos (ndj) = 0 9j = 


2j — In 


(2.4) 


n 


for j = 1,2, 


, n 
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2.3. Extrema 

|T n (x)| = 1 is the extrema of T n (x ) 


\T n (Vk)\ = 1 7? fe = cos^— J for A: = l,2,---,n 
T„(±l) = (±1)” 


(2.5) 

( 2 . 6 ) 


2.4. Polynomial Expression 

A simple method for finding Chebyshev polynomial expressions can be done using 
trigonometry and a recurrence relation. Start with the first two polynomials, 
T 0 (x) = cosO = 1, and T\{x) = cos 9 = x. Prom the relation 

2 cos 9 cos nd = cos (n + 1) 9 + cos (n — 1) 9 (2.7) 

2 xT n {x) = T n +\[x) + T n -i(x) 

T n+ 1 {x) = 2xT n (x) — T„_i(x) (2.8) 

T 2 (x) = 2xTi (x) — T 0 (x) = 2x 2 — 1 

T 3 (x) = 2 xT 2 (x) - Ti(x) = 4x 3 - 3x 

This method of arriving at the n th degree polynomial is preferred when writing 
computer programs [12]. 

The following is a more rigorous method for determining Chebyshev Polyno- 
mials of any degree. First write out DeMoivre’s Theorem. 

e 1 9 = cos 9 + i sin 9 

e ind = (cos 9 4- i sin 9) n — cos(n0) + i sin (n9) (2.9) 


by the binomial expansion, 

(cos 9 + i sin 9) n = cos” 9 + T" j cos” -1 9{i sin 9) + 


71 cos” -2 * 


(isin0) 2 -| 1" (n)(* S in^)” ^-10) 


Equating the real parts of (2.10) gives 

cos n9 = cos” 9 -1- cos” -2 9(i sin 9) 2 + cos” -4 9(i sin 9) A 


-\ + cos" 2v 9{i sin 0) 2v 


( 2 . 11 ) 
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cos n6 — y : Aj y: B kj — A? [B 0j + By H h By] 

j =0 fc =0 j =0 

= Aq [Boo] 

+Ai [Box + Bn] 

+A2 [B02 + B12 + B22] 

+Aq [B 0? + Big + • • • + Bgg] 

+Aj [Bog 4- B2 V + ■ • • + Bffrj] (2.15) 
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Factor the Bkj terms 


Boo [Ao] + Boi^ 4 l + B02A2 + • • * + BorjArj 
+B\\A\ + B12A2 + • • • + Bir,A v 


+B qq A q + B qq+ iA q+ i + B qq+ 2A q+ 2 H I- BqqA-q 


if]— 1 Ar}— i + Brj—irjAjj 

H~ BriqA-q 
V 

= Y, BjkAk for j = 0, 1,2, — ,77 

k=j 

77 7 ) 

- £ £ B ikA k 

3 - 0 k=j 

Substituting back in the expressions for Ak&nd Bjk gives 


v *7 

EE 

j =0 k—j 


(-i y ( . j cos 2j 0 


k-D A 


cos 


n—2k , 



+ (-1)" { n cos n - 2n e 

v 2r ?/ \3j 


n 


-E ( - 1)2m ( 2 (. + ? ) A j J 


W + A 


COS 


~ 2q 9 for q = 0,1, 


, n 


(2.16) 

(2.17) 


(2.18) 
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showing that (— 1) 2 -* +9 = (— l) 2j (— l) 9 = (— l) 9 gives the final expression for 
the Chebyshev polynomial of degree n. 

T n (x) = cos n6 = 5I(-l) fe 2 j ) (fc) cos n_2fc e (2.19) 

substituting for x = cos 6 


T n (x) 

t n -(2k+l) 


tn— 2k 


to + t\X + • ■ • + t n x n 
0 


(-D*£ 


j=k 




for k = 0, 


n — 1' 
2 . 


for k = 0, • • • 


n 

2 


( 2 . 20 ) 


The above equation states that if n is even then all powers of x are also even 
and if n is odd then all powers of x are odd. The first few polynomials are given 
below. 


To(x) = 1 T 2 (x) = 2x 2 — 1 T 4 (x) = 8x 4 — 8x 2 + 1 

T\(x) = x T^(x) = 4x 3 — 3x T 5 (x) = 16x 5 — 20x 3 + 5x 

2.5. Methods of Summation 

In the next section, certain recurrence relationships involving derivatives of Cheby- 
shev polynomials will arise. These recurrence formulas will give a rather compli- 
cated summation expression. To properly deal with these summations, a technique 
know as methods of the summation [11] will now be introduced. 

To begin, we first define the difference operator 

A/ (x) = / (x + h)-f (x) (2.21) 

and let h = 1. 

A/(x) = /(x + l)-/(x) 

The difference operator A is a linear operator and behaves 
exactly) to a differential operator. 

A [ af (x) + bg (x)] = aA/ (x) + bAg (x) (2.23) 


( 2 . 22 ) 

similarly (but not 
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A[/(x)s(x)] = /(x + 1)j(x+1)-/(x)j(x) 

-[/(x+l) S (x)-/(x + l) S (x)] (2.24) 

= /(x + l)(j(x + 1) -g(x)] + j(x)[/(x+ 1) -/(x)] 

A [/ (x) g (x)] = / (x + 1) Ag (x) + g (x) A/ (x) (2.25) 

Similarly 

a /(s) = 9 (£) A/ (x) - / (x) Ag (x) ( 2 _ 26 ) 

ff(i) 3^ + l)</W 

Here are some examples. 

Ac = 0 
Ax = 1 
Ax 2 = 2x + 1 
Ax 3 = 3x 2 + 3x + 1 


where c =constant and n is a positive integer. 

Repeated operations give 

A 2 x 2 = A (Ax 2 ) = A (2x + 1) = 2 
A 3 x 2 = A (2) = 0 

in general A n cx n = cn\ and A m cx m = 0 for m > n. 

Now we shall introduce the anti-difference operator A -1 

if A g (x) = / (x) 


then 


g( x ) 

AA ~ l f{x) 

f( x ) 


A-V(x) 
A g ( x ) 

A g (x) 


Note that AA -1 = 1 but A -1 A ± 1! That is to say the operations do not 
commute. 


A V(x)] = / (as) 


where c is some constant. 


A 1 [A/ (as)] = / (x) + c 


(2.27) 
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Example 2.1. Verify that A _1 x = | (x 2 — x) + c 


AA _1 x 



Now we will introduce some notation on factorial products. 

X (n) = x (x - 1) (x - 2) • • • (x - n + 2) (x - n + 1) (2.28) 

x^") is called the n th factorial of x. Below is a list of the first few factorials. 


x® = x (x — 1) = x 2 — x 

x (3) = x (x - 2) (x - 1) = X 3 - 3x 2 + 2x 

the difference operator of these factorials are 

Ax^ = Ax = 1 

Ax^ _ ^ ^j.2 _ — 2x = 2x^ 

Ax ( 3) = A (x 3 - 3x 2 + 2x) = 3 (x 2 - x) = 3x (2) 
in general Ax^ = 

The n th factorial of x, x^ n K can be given in terms of a polynomial expression. 
I (n) = V" 2 + - • •+-s< 1) x+5f = jr Sl m) xW (2.29) 

771=0 

the set S'4 1 \ • • • , 5^} is known as Sterling’s Numbers of the first kind. 

Likewise, any power of x can be represented by a polynomial of factorials. 

x n = e^x^ + ef-W— 1 ) + ■■■ + + ei 0) = (2.30) 

771—0 

The set 1©^, ©£>, • • • , 6^ n) } is known as Sterling’s Numbers of the second kind. 
An excellent source for fisting a large number of Sterling’s Numbers is given in 
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the Handbook of Mathematical Functions [6]. A list of the first few of Sterling’s 


Numbers of the first and second kind are 

given below. 

x<°> = 1 

II 

T— t 

II 

o 

X^ = X 

x 1 — X^ 

x< 2 > = x 2 — X 

x 2 = x< 2 > + x^ 

x( 3 > = x 3 — 3x 2 + 2x 

x 3 = x^ + 3x^ + x^ 

x (4) = X 4 - 6x 3 + llx 2 - 6x 

x 4 = x^ + 6x^ + 7x^ + x 

x (5) = X 5 - 10x 4 + 35x 3 - 50x 2 + 24x 

x 5 = x (5) + 10x (4) + 25x (3) + 15x (2) + 


The introduction of the difference operator and the factorial notation will now 
be used together to show a convenient way of solving for the limits of summations. 
Look at a simple combination of the difference operator and a summation sign. 

£ a/ (*>=£[/ (*+ i)-/w] 

2—1 2=1 

= I/(5) + /(4) + /(3) + /(2)]-[/(4) + /(3) + /(2) + /(l)] 

= / (5) -/(!)= /Ml”? 


in general 

£ A/ (x) = / (x) 

x—a 

Now let A / (x) = g (x) 

/ (x) = A ~'g (x) + c 


x=6+l 


(2.31) 


b b 

2=6+1 


2=6+1 

Ea/W = Y,g(x) = f(x) 

= A 1 g(x) + c 

II 

> 

1 

<cT 

x =a x—a 

Therefore 

b 

x =a 

2=6+1 

x=a 

Es(x) = A ^ g (x) 

x=a 

x =a 



2 = 64-1 


(2.32) 


This formula gives a method of taking a function under summation and finding 
its limit assuming one exists and its anti-difference operator can be found. 
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Example 2.2. 


5> 2 = £(x< 2 )+xW) 

X=1 1=1 

x( m ) = A 1 ^* 1 so = Aj and x^ = A^- 

x=n-f 1 
x=l 

_ (n + l) (3) (n + 1) (2) I® 1( 2 ) 

3 + 2 _ 3 ~ 2 

(n + l) (3) (n + 1) (2) 

3 + 2 

(n + 1) (n) (n — 1) (n + 1) (n) 

_ - + - 

^ 2 _ (2 n + 1) (n + 1) n 

h x " 6 

Section 2.1 through 2.4 introduced the definition and some of the properties 
of Chebyshev Polynomials as well as a derivation of the polynomial itself. Sec- 
tion 2.5 showed a nice trick for evaluating summations which have limits and an 
antidifference operator. The definition and properties of Chebyshev Polynomi- 
als where given to aid anyone interested in more detailed derivations. Section 
2.5 will be used in section 3, particularly section 3.2, where recurrence relation- 
ships between functions expanded in terms of Chebyshev Polynomials and linear 
differential operators of these functions, will be dealt with. 

2.6. Problems 

Exercise 2.1. Find the first four powers of x in terms ofT n (x ) by a recurrence 
relation. 

Exercise 2.2. Prove that T n (x ) is commutable; that is T r (T s (x)) = T rs (x) = 
T st (x) . Note that there exists one and only one polynomial which is commutable 
with T n (x). Therefore if a polynomial p n (x) is commutable with T n (x) then 
Pn(x) = T n (x). 
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Exercise 2.3. Prove that the coefficient on the x n term of the Chebyshev poly- 
nomial is 2 n ~ 1 . Hint:( 1 + x) n = £” =0 (”)^ for n > 0, n a positive integer. 


Exercise 2.4. Show that £" =1 


n 2 (n + l) 2 
4 


3. Differentiation of Chebyshev Polynomials 

This section will illustrate some of the derivations and applications of the deriva- 
tives of Chebyshev polynomials. One of the most important applications is in 
approximating a function. We will then proceed to take derivatives of the func- 
tions and show recurrence relationships between the coefficients of the functions 
and its derivatives. This last step is important in the application of spectral 
methods. 


3.1. Stiirm-Liouville solution and Orthogonality 


Chebyshev polynomials can also be defined as the function which satisfies the 
following differential equation. 


(i 



d 2 y 

dx 2 


- x^j- + n 2 y = 0 

dx 


(3.1) 


Since (3.1) is a Stiirm-Liouville problem, the Chebyshev polynomial is a par- 
ticular hypergeometric function. 

Chebyshev polynomials are orthogonal functions with respect to a weighting 
function of (1 — x 2 ) 2 . 


j l T n (x)T m {x ) (l- 

where 



^ dx — ^ Cn^nm 

(3.2) 

n — 0 
n > 0 

(3.3) 

n^m 
n = m 

(3.4) 


Orthogonality in combination with its geometric rate of convergence are the rea- 
sons why this function is so useful. 
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3.2. Recurrence Relations 

Let’s say we wanted to approximate a function by a finite series of Chebyshev 
polynomials. 

N 

u (x) ~ U N ( X ) = 53 OnT n {x) (3.5) 

n= 0 

say we also wanted to approximate the derivatives of the function. 

U N ( X ) = 53 a n )T n( x ) (3-6) 

n= 0 

The derivative of u is related to N Chebyshev polynomials by a separate set of 
coefficients K By taking the derivative of (3.5) we obtain. 


N 


U N ( X ) — 53 a ‘ n 'I'n{ x ) 


1 


(3.7) 


Note that now the summation begins with n = 1 and not n = 0, while the 
summation of (3.6) starts at n = 0. 

A relationship between a n and can be found. To prove this, we will make 
use of the following trigonometric identity. 


2 sin r# cos p# = sin (p + r) 6 — sin (p — r) 9 
sin (p + r) sin (p — r)6 


2 cos pd = 


sin r0 


sin r9 


Now let p = n and r = 1 


2 cos nO = 


sin (n + 1) 9 sin (n — 1) 6 


sin 9 sin 6 

which upon substitution of the definition of T n (x ) and T' n {x) gives. 

2 T n (x) = ^±iM _ Zk r . jfr) 

V ' n + 1 n- 1 

Substituting (3.9) into (3.6) and equating it to (3.7) yields. 


N - 1 

5> 

n= 0 


w± 


n+1 n — 1 


= 53 a *X(x) 


n= 1 


(3.8) 


(3.9) 


(3.10) 
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With the knowledge that T_ m (x) = T m ix), we can start writing out terms of 
(3.10). 


a, 


(i) 


T[(x) Hit r) 


-1 


+ 


,n fTl,(x) TUlM 

2^ a n 


n=l 
N 


n + 1 
= 2 a n T n {x) 


n — 1 


71=1 


2 aQ^Ti(x) + a[ 


(i) 




r,(x) _ r n (x) 
2 0 


+«?> 


ZkM _ r i(*> 

3 1 


ri(») _ T,(x) 

4 2 


+ • • • — 2a 1 T 1 + 02-^2 (^O ~b 2a 3 r 3 (ar) + 


(3-11) 


(3.12) 


Since { T n (x ) : n = 1, 2, • • • , N} make up the basis of an N dimensional Chebyshev 
space, the coefficients on each term can be set equal to each other. Also note that 


77i — ►() 772 m-+0 

sin mO 0 

sin 6 sin 6 

(3.13) 

2 ag' - a? 

— 

2a\ 


4 1} 4 1} 
2 2 

= 

2d2 


4 1} 4 :) 

3 3 

= 

2az 


a W - a W 
a k a k + 2 

= 

2 (fc H- 1) Ufc+i 

(3.14) 

07 (1) (!) 

for 1 < A: < N 

(3.15) 


in general 


where Ck is defined in (3.3) 

The same relationship can be developed for the ( q — 1) £A derivative of (3.9). 




T“ x (x) T™,( i) 


(9) 


n > q 


(3.16) 


n + 1 n — 1 

where T^{x) is the q th derivative of T n (x). Applying the same procedure to (3.16) 
would yield 

cJp-a^ = 2(k + l)a^P 0 <k<N-q (3.17) 
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Starting with (3.14) a simplified form of the recurrence relation can be devel- 
oped. Notice that a* = 0 for k > N. Start with k = N — 1 and work in descending 
order. 

2 (N) a^i 

2 ( N — 1) ajv_:L 

2 ( N — 2) a^v-2 

c-id^ — = 2(3 )a 3 

Ci<2^ — = 2(2 )a 2 

Coa^ — = 2(1 )a 2 (3.18) 

Simphfying and substituting yields 

fljv-i = 2 (N) aw 

a !v _ 2 = 2 (JV — 1) a^_i 

a$_ 3 -2{N)a N = 2{N -2)a N _ 2 

®iv-4 — 2 (iV — 1) a^v-i = 2 (iV — 3) ajv-3 

dj ^ = 2 • 2a 2 -I - 2 ■ 404 + • ■ ■ + 2 (iV — 1) <i/v— 1 

coOq ^ = 2 • 2al + 2 • 4<i3 -(-••■ + 2 (iV) (3.19) 

Put (3.19) into summation form. 

N 

c n ag ) **2 £ P a p ( 3 -20) 

p=n-f 1 
p+n odd 

The relationship between the q th and the q th — 1 derivative coefficients can also 
be written. In addition (3.20) does not need to be truncated to N terms. The 
combination of the two gives. 

OO 

= 2 P a p~ l) (3-21) 

p=n+l 
p+n odd 
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To carry this exercise further, let’s say we were interested (which we are) in 
the second derivative recurrence relation between a n and affi. Let q = 2 in (3.21) 
and substitute (3.20) in for a^. 


Cjid. 


i? = 2 £ 


m 


m=n- j-1 
m+n odd 


cy OC 

— E P a p 

^ rn p=m -\- 1 
p+m odd 


(3.22) 


Since m will always be greater than zero, Cm = 1. Simplifying (3.22) gives 


CnO 


(2) 


oo 

= 4 E m 

m=n- j-1 
m+n odd 


oo 

E P a p 

p=m - hi 
j)+m odd 


(3.23) 


To get rid of the oo in the inner summation we can switch the order of sum- 
mation. Note that the summations can only be switched when the summations 
are uniformly convergent. Start by expanding the inner summation. 


OO 

Cnc4 2) = 4 m[(m-|- 1)^+! + (m + 3)a m+3 -l ] (3.24) 

m=n- f-1 
m-bn odd 


Followed by the expansion of the outer summation. 

icna^ = (n -I- 1) [(n + 2) a n + 2 + (n + 4) a n+ 4 + (n + 6) a n + 6 + • • •] 
+ (n + 3) [(n + 4) a n + 4 + (n + 6 ) a n +6 + • • •] 

+ (n - 1-5) [(n + 6 ) a n+ 6 + (n + 8 ) a n+ 8 H ] 


-i- (n + 2i — 1 ) [(n + 2 i) a n+ 2 i + (n + 2 i + 2 ) a n+ 2^+2 + • • '13.25) 

Regroup the a t terms and their corresponding coefficients. 

4 {(n -I- 2) a n+2 [n + 1] 

+ (n + 4) a n+ 4 [(n + 1 ) + (n + 3)] 

+ (n + 6 ) a n+6 [( n + l) + {n + 3) + (n + 5)] 

+ (n + 2 i) a n+ 2 i [(n + 1 ) + (n + 3) 4 h (n + 2i — 1 )]} 

: (3.26) 
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Write the inside brackets as a summation. 

n+2i — 1 

4 (fl + 2i) Q“n+2i £ m for t = 1,2,3, ••• (3.27) 

m=n+l 
771 -f-n odd 

and finally 

OG p — 1 

CnQ^ 2 - ’ = 4 Y2 P a P Y1 m ( 3 - 28 ) 

p=n +2 m=n+l 

p+n even m+n odd 

Recalling from section 2.5 on the methods of summation, we saw ways of 
finding the value of 5Z P mI n +i m ■ ^ would be more aesthetic and computational 

m+n odd 

less intensive if there existed a single value that would replace this summation. 
Begin by expanding the summation. 


p-i 

y. m = (n + 1) + (n + 3) H h(p— 1) (3.29) 

m=n+l 
m+n odd 


2-n 

2 


E 2 — 


? 2 


= Y, n + (2® - 1) = y (n - 1) + 2 y x 


(3.30) 


2=1 


2=1 


2=1 


Let k = ++ . Look at first of the two summations on the right hand side of 
(3.30). 

k k 2=fc+l 

y (n — 1) = (n — 1) = (n — 1) x 

2=1 2=1 X =1 

Now look at the second summation of (3.30). 


= (n — 1) k 


( 2 ) 


25> = 2^aV = 


= + 2 > 


2=1 


2=1 


2=fc+l 


= {k + l)k 


i=i 


p-i 


y m = (n — 1) k + (k + 1) k 


771 = 71+1 

m+n odd 


= k [(n — 1) + (k + 1)] = k (n + k) 

f p — n\ ( p-n\ 1 ( 2 2 \ 

—) { n+ ~2-) = -A p ~ n ) 


(3.31) 


(3.32) 

(3.33) 
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(3.34) 


Substituting back into (3.28) and truncating to N terms gives. 

Cnd { n ’ = Y P(P 2 ~ n2 ) a P 

P=Tl+ 2 
p-f n even 

Equation (3.34) is a convenient way of representing the second derivative of a 
function in terms of a series of Chebyshev polynomials. 

There is another formula which is very important in solving boundary value 
problems. The equation gives the value of the p th order derivative of a Chebyshev 
polynomial when evaluated at one of its end points, 1 or — 1 . 

^r„ (± D = (±ir* n ■ <3 ' 35) 

Equation (3.35) is given without proof but has been verified by the author in 
several problems. Notice that instead of using coefficients to relate the derivative 
of the function to the zeroth order derivative of the Chebyshev polynomial, one 
directly takes the derivative of the Chebyshev polynomial and then evaluates it at 
the end points. Equation (3.35) gives an extremely useful method for evaluating 
derivatives in boundary conditions. 

3.3. Matrix Notation 

In the last section, summation notation was used to develop recurrence relation- 
ships between expansion coefficients of a function (i.e. a n ) and the expansion 
coefficients of derivatives of the function (i.e. a^). In this section, a matrix nota- 
tion will be used to develop the same relationship between a n and as well as 
other expansion coefficients. In my personal experience this method proved the 
most insightful and easiest to derive. It does not, however, deal very well with 
recurrence relationships such as the following expansion. 

OC 

xu N (x) = Yb n T n (x) (3.36) 

n= 0 

or 

xu N (x) = Yb ( n ) Tn(x) (3.37) 

n=0 

The matrix notation will be used in the next section which explores the application 
of Chebyshev Polynomials to spectral methods. 
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Look back at the recurrence equation (3.18). Put the equations into matrix 

form by letting a = (a^, . . . , a $_ : ) and remembering that cq = 2, Cn = 1 

for n > 0. 


( 0 

0 

0 

U 

which can be 


Since C is an upper triangular matrix, det (C) is the product of its diagonal entries 
which here det (C) = 2. As det (C) ^ 0, C is non-singular and C -1 exists. This 
allows (3.38) to be multiplied by C -1 . 

a (1) = C -1 A a (3.40) 

Let C -1 • A = E . This gives a^ 1 ^ = E • a. In general the expression for the q th 
derivative can be written 

a (?) = E-a (9 - 1} (3.41) 

a (2) = E • a (1) = E E a = E 2 -a (3.42) 

in general we have 

a w = E fc -a (3.43) 

where E fc =E E E E 

n ' 

k times 


( 2 0 
0 1 
0 0 


VO •• 
2 0 
0 4 

0 0 


0 

0 


-1 0 ••• 0 \ 

0 -1 ••• 0 

1 0 0 

1 0 

o 

... o\ 

0 

0 

2(iV - 1) 

0 2 N ) 


written in matrix form as 

C • a^ = A • a 


\ 


(i) 


a 


\ 4-i ) 

l a 0 ' 

a i 


&N - 1 
V a N ) 


(3.38) 


(3.39) 
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It can be shown that 

0 1 0 

0 0 4 

0 0 0 


0 5 0 7 0 ••• 2 N 

8 0 12 0 0 0 

0 10 0 14 0 ■■■ 2N 


(3.44) 


E 2 = 


'n+1) ( m 

+i)j - 

« 1 
II 

m 

for 

n = 
m - 

= 0,1, 
= n + 

■ ■ • N — 2 
1, n + 3, • • • , N 

(3.45) 

0 

0 

4 

0 

32 

0 

108 

0 

256 \ 


0 

0 

0 

24 

0 

120 

0 

336 

0 i 


0 

0 

0 

0 

48 

0 

192 

0 

480 ; 







80 

0 

280 

0 i 

(3.46) 







120 

0 

384 • . : 



(3.47) 


which corresponds with 

E = [E ( „ +1) (*.,)] = fp (p 2 - O n p Z “ + 2,’„ + 4 . 7. 3 J V < 3 ' 47 > 

The main focus of this section was to take a function and derivatives of that 
function, which we wanted to expand into Chebyshev Polynomials, and show 
that the expansion coefficients of each can be related. This was done by two 
methods, summation notation and matrix notation. The matrix notation is useful 
in demonstrating how to solve for systems of differential equations. It will also be 
used to solve for a variety of eigenvalue problems. The recurrence relationships 
are preferred when only the solution vector to a problem is desired. Most of 
the common recurrence relationships for linear operators are tabulated in the 
appendix of Gottlieb and Orszag’s book, Numerical Analysis of Spectral Methods 
[41. 
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3.4. Problems 


Exercise 3.1. Prove that T n (x) = cos nO is the solution to 


i 1 ~ x 2 ) 

Exercise 3.2. Prove that 


d 2 T n (x ) dT n (x) 


dx 2 


x- 


dx 


+ n 2 T n (x ) = 0 


J ^ T n (x)T m {x) (l - x 2 ) 2 dx = ^CnS n 


Exercise 3.3. If <p r (x) is a real orthogonal function with respect to a weighting 
function w (x) and w (x) > 0 on the interval X = —1 < x < 1. Prove that 0 r (x) 
has r real and distinct roots on the interval X. 


Exercise 3.4. Prove 


1 N 

<v4 3) = 7 E - l] [(m + nf - l] a m 


(3.48) 


771 = 71+3 
m+n odd 


Exercise 3.5. Prove 


N 


<v4 4| = ^7 E P P 2 (p 2 -4) -3nV + 3nV-n 2 (n 2 -4)' 


a, (3.49) 


p=n+ 4 
p+7i even 


Exercise 3.6. Prove that 
its interval is 


the derivative of T n (x) evaluated at the boundary of 
£T„(±l)=(-l) n+ ‘n 2 


4. Spectral Tau Methods 

4.1. Introduction 

Spectral Methods axe a particular numerical scheme for solving differential equa- 
tions. It is a discretization scheme developed from the Method of Weighted Resid- 
uals (MWR) [1], The Tau method is one of the three most commonly used tech- 
niques of spectral methods. The three techniques are Galerkin, collocation and 
tau. This section will concentrate solely on the tau method. The interested reader 
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is referred to Spectral Methods in Fluid Dynamics [2] for examples of the other 
two methods. 

Before further explanation of the spectral tau method, a brief review of the 
method of weighted residuals may be appropriate. Given the problem 

^ = L u + f (4.1) 

at 

Bu = 0 (4.2) 


where L is some linear spatial operator and B is a linear boundary operator. We 
want to express u (x, t) as an infinite sum of some trial function 4> n (rr). 

OO 

u(x,t) = ^2a n (t)4> n ( x) (4.3) 

n= 1 

the function is then approximated by truncating the series to N terms. 

N 

u N (x,t) = ^2a n (t)<t> n (x) (4.4) 

n— 1 

Now choose a test function ip m (a:) such that 4> n (x) is orthonormal to ib m ( x ) on 
some interval 2 — a < x < b for a given inner product (*,*). 

rb 

(<Pn (^/) j 1pm (®)) = / (•*-) 1pm (^-) dx b nm (4-5) 

Ja 

where ip m (x) is the complex conjugate of ip m (x)- The essence of the Galerkin 
method is to find the set {a n (t) : n = 1, 2, • • • , N} such that the errors {e n } are 
minimized. The error is defined as 


du 

dt 


du n 


Lit + Lu n f -b fn — e n 


To find the set of coefficients {a„ (t) : n — 1, 2, • • • , N} multiply both sides of the 
governing equation (4.1) by the test function and take the inner product {*,*) 
such that 


al (VVn 0*0 > U N ) — (ipm (%) > Lujv) + (ipm ( x ) > /) (4-6) 

at 

N rfr, — N 

+ (4.7) 

n=l at n = 1 
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Take the complex conjugate of the equation and note that (ipm, <Pn) = <5 nm . This 
yields 

(*) , L0n (x)) + {/, iftm (x)) (4.8) 

ai n= 1 

Equation (4.8) can be solved by well known methods which will depend upon the 
operator L and the inner product. 

4.2. Spectral Tau Methods 

Notice that the choice of the trial function 4> n (x) must satisfy the essential bound- 
ary conditions. For particular problems, the boundary conditions can be compli- 
cated. Choosing a trial function which satisfies the boundary conditions exactly 
can become extremely difficult if not impossible. The spectral tau method offers 
a solution to this problem. The solution is quite simple in its concept. Instead 
of trying to satisfy the boundary condition by choosing the proper trial function, 
just add more constants to the expansion. 

N+k 

u N (x, t) = a n (t) <j) n (x) (4.9) 

n=l 

Here N is the dimension in the domain and k is the number of boundary condi- 
tions. So before where <p n (x) was forced to satisfy the boundary conditions, we 
just add k more coefficients, {a n (t) : n = N + 1, • • - , N + k}, to be solved. The k 
additional equations needed are produced from the boundary conditions. 

N+k 

53 a r^4>n (x) = 0 (4.10) 

n=l 

So the complete set of N+k equations to be solved for the set of N+k coefficients 
are the N equations of (4.8) and the k equations of (4.10). 

The application of Chebyshev polynomials to the spectral tau method is to let 

4> n (i) = T n (x) (4.11) 

and 

^m(x) = — T m (x) (l -x 2 ) 2 (4.12) 

■nCm v ' 

The best way to explain the use of Chebyshev spectral tau method is to work out 

an example. 
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Example 4.1. The following system is a second order differential equation with 
homogeneous boundary conditions. 


eui a. 

to| £ 

II 

(4.13) 

u (— 1) = 0 

(4.14) 

du( 1) 

dx 

(4.15) 

Solve for u using the Chebyshev spectral tau method. This system has the solution 

u (x) = \x 2 — x — |. 


First apply 

N 


U ( X ) = J2 a nTn{x) 

(4.16) 

n~ 0 


Here the order of the summation has been changed to 

n = 0 to n = N. We will 

need the first and second order derivatives of u (x) as 
this, equations (3.20) and (3.34) will be used. 

they relate to T n {x). For 


(4.17) 


(4.18) 


Substitute (4.16) into the governing equation (4.13) and the boundary conditions 
(4.14) and (4.15). 


£>r„M = i 

71=0 

N 

(4.19) 

£a n T n (-l) = 0 

71=0 
N- 1 

(4.20) 

E 4 ,: >r„(i) = o 

71=0 

(4.21) 


Alternatively and preferentially, we can use equation (3.35) for (4.21). Equation 
(3.35) is preferred because it is easier to program in FORTRAN. When (3.35) 
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is used on the boundaries, the lower index must be raised to the order of the 
derivative. So in general, 


d p u N _ ^ dPT^jx) 
dxP i^ p am dxP 


(4.22) 


For this problem we will use the recurrence relationship of a), 1 ). In the next and 
subsequent problems, (3.35) will be used. As an exercise, the reader should use 
both methods to verify that each will obtain the same result. So (4.21) can be 
written as 

E a nT' n (1) = J2 ^ ( 4 - 23 ) 

n=l n= 1 

The first step will be to take the inner product of the domain equation (4.19) 
E °n 2) (Tn(x),1pm{x)) = (l,V>m(*)> = (T 0 (x) , l/> m (x)) (4.24) 

n—Q 

where 

(T n (x), ip m {x)) = — [ T n (x)T m (x) (l - x 2 ) 2 dx = 6 nm (4.25) 

7 TCn J - 1 v 7 

this gives. 

(2) c f 1 n = 0 
°n -<5on = | Q Q<n < ^.2 

On the boundaries, we need only evaluate the Chebyshev polynomials, T n (l) = 1 
and T n (— 1) = (— l) n . Upon substitution into (4.20) and (4.23) yields 

EM- l) n = 0 (4-27) 

71 = 0 

N 

E n 2 a n = 0 (4.28) 

n= 1 

This gives N + 1 equations, N — 1 for (4.26) plus the two boundary equations 
(4.27) and (4.28). Let N = 5 and write out equation (4.26) using the recurrence 
relationship for the second derivative expansion coefficient from equation (3.34). 

<2 q ^ = — [2 • 4ci2 + 4 • 16(24] = 1 
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(4.29) 


<zf = 3 -(9 -l)o 3 + 5 -(25-1 )o 5 = 0 
af = 4 -(16-4)a 4 = 0 
of - 5 -(25 -9) o 5 = 0 

The boundary conditions are 

00 — 01 + 02 — 03 + 04 — 05 = 0 (4.30) 

Oi + 4 o 2 + 903 + 1604 + 25a 5 0 (4.31) 

Put the N + 2 equations into matrix form. 

/ 0 0 4 0 32 0\/oo 

0 0 0 24 0 120 Oi 

0 0 0 0 48 0 a 2 

0 0 0 0 0 80 o 3 

1- 11-1 1 -1 04 

0 1 4 9 16 25 / \ o 5 

The solution to (4.32) is a = (— |, -1, 1,0,0, o) T , substitute into (4.16). 

u N (x) = -^To(x) -Ty(x) + ^T 2 (x) 

= — - — x + - (2x 2 — l) = \x 2 — x — ^ 

4 4 1 ' 2 2 

which is the exact solution. 

The matrix notation is a heuristic lesson in solving differential equations with 
the spectral method. The matrix notation is also an excellent way of solving for 
eigenvalue problems. When only the solution vector is needed, the use of the 
recurrence relation (3.17) is preferred. The reason is purely computational. In 
general, solving eigenvalue problems take of order N 3 operations. The recurrence 
relation will usually form a tridiagonal system (for one dimensional problems) 
which can be solved with much more efficient algorithms [12]. An example of how 
this is done is given in Spectral Methods in Fluid Dynamics pp. 129-131 [ 2 ]. 

Now we will demonstrate the use of Chebyshev spectral tau methods in solving 
eigenvalue problems. 
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Example 4.2. This example deals with an eigenvalue problem with homogeneous 
boundary conditions. We will try to solve for the first few eigenvalues, A. 


d 2 u 

dy 2 


+ \ 2 u = 


U\ 


y=o 


du 

dy 


+ u 


v=l 


0 

0 


= 0 


ly=l 


The analytical solution for A is 


(4.33) 

(4.34) 

(4.35) 


A = — tan A (4.36) 

In this example, we will more closely look at the details of each step. This 
careful analysis will shed light on some of the intricacies involved in this method. 
There are three key steps that must be applied between the system of equations, 
which comprise of the governing equations and the boundary conditions, and the 
set of equations that need to be solved for the expansion coefficients, a n . The first 
step is to interpolate or approximate the function u (x) into Chebyshev space. 


u 


N 

(x) — un (x) = ^ ^ a n T n {x) 


n—0 


(4,37) 


The interpolated function, u N (x), is projected into the equation space. This step 
will be explained in more detail later. And finally, take the inner product of the 
domain equations (governing equation) and evaluate the projected functions on 
the boundaries. 

Notice that the Chebyshev polynomials he in the interval of —1 < x < 1 
and the problem is stated on the interval of 0 < y < 1 . In order to transform 
the problem into the proper space, we will apply a stretching function (more 
commonly called a coordinate transformation) . The stretching function we need 
here is x = 2y — 1 . Transforming the system of equations (4.33) through (4.35) 
yields. 


' dx\ (fu , 

dy) dP + U 


u 


x=-l 


1 dx\ 

du 

+ u 

\dy 7 

dx 

X=1 


0 

0 


(4.38) 

(4.39) 

(4.40) 
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X=1 


0 



(4.41) 

(4.42) 

(4.43) 


Using the derivative — = 2 and interpolating u (x) into (x) will give. 

dy 


d?UN 
4 dx 2 


4- A 2 ujv 


U N | I= _l 


2 


dux 


dx 


+ Upf 


x=l 


X— 1 


0 

0 

0 


Notice that now the problem is cast into u N (x), the solution technique will find 
the exact solution to u N (x). If u (x) is exactly cast into u N (x) then the exact 
solution to u (x) will be found. In general if u (x) is a polynomial of degree equal 
to or less than the degree of the Chebyshev polynomial being used, then the exact 
solution to u (x) will be found. That is why example 4.1 came up with the exact 
solution. In this problem, and in general, we will not find the exact solution, since 
there exists an infinite number of eigenvalues. 

Substitute (4.37) for ti#(x) into (4.41), (4.42) and (4.43). 

4xfa^T n (x) + A 2 f^a n T n (x) = 0 

rt=0 n=0 

£> n r n (- 1) = o 

71=0 

2f> 2 a n T n (l) + £xT n (l) - 0 

n=l n=0 

The last key step is to take the inner product of (4.44) and evaluate T n (— 1) and 
T n ( 1) in (4.45) and(4.46), respectively. Remember from the Tau method, that 
the inner product of (4.44) will be taken N - 2 times; N minus the number of 
boundary conditions. It is important to note that we will subtract from equations 
in which the boundary conditions apply to. This will become important when 
there is more than one differential equation. 


(4.44) 

(4.45) 

(4.46) 


4a^ 2) 4- A 2 a n = 0 for n = 0, 1, ■ • • , N — 2 (4.47) 

f> n (-l) n = 0 (4-48) 

71=0 

2 £> 2 a„ + £>„ = 0 (4.49) 

n= 1 n=0 
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By using the matrix notation developed in the last part of section 3, (4.47) through 
(4.49) can be simplified. 

(4E 2 + A 2 l)-a = 0 (4.50) 

(— 1)" / • a = 0 (4.51) 

(2 n 2 / 4- /) • a = 0 (4.52) 

where (-l) n /= (l, -1, 1, -1, . . . , (-1)*) and 2n 2 I= (0,2, 8, 18, . . . ,2N 2 ). Note 
that 4E 2 + A 2 I is an (TV — 1) x (N + 1) matrix. The last two rows needed to make 
a square matrix come from (4.51) and (4.52). Also notice that the I in (4.51) is a 
vector with N + 1 entries equal to one. This notation leads into the format that 
will be needed to write a numerical computation. 

Now the initial statement of the problem was to find the first few eigenvalues 
of the system. We can achieve this by rewriting (4.50) through (4.52) as. 

4E 2 ■ a = — A 2 I • a 
(— l) n I ■ a = 0 
(2 n 2 I + /) • a = 0 

Equations (4.53) through (4.55) give a generalized eigenvalue problem 

A • a = uC • a (4.56) 

where for N = 5 

0 0 16 0 128 0 

0 0 0 96 0 480 

0 0 0 0 192 0 

0 0 0 0 0 320 

1-11-1 1 -1 

1 3 9 19 33 51 



4E 2 \ 

(- 1 )”/ = 
2 n 2 I + / / 


(4.53) 

(4.54) 

(4.55) 



/ — 1 0 0 0 0 0 \ 

0-10 000 
0 0 -1 0 0 0 

0 0 0-100 
0 0 0 0 0 0 

0 0 0 000 / 


(4.58) 
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and u — A 2 is the eigenvalue. Note that both A and C are (N + 1) x ( N + 1) 
matrices. If the matrices A and C were imported into a computer program and 
then subsequently passed to a generalized eigenvalue solver, the first few values 
of A can be found. 

To calculate the first few eigenvalues, I used MathCad version 5.0+ . A 
general eigenvalue solver is available in this software. Below is a table of the first 
four eigenvalues for different values of N and the analytic values of A found using 
equation (4.36). 

Table 4.1 First four eigenvalues of Example 4.2 
for various values of N and the actual value. 


N 

4 

5 

6 

7 

8 

9 

Actual 

A: 

^2 

^3 

a 4 

The next exa 

1.990 

5.73 

imple v 

2.029 
4.966 
9.377 
22.44 
411 illus 

2.029 
4.917 
8.204 
13.96 
trate h< 

2.029 
4.914 
8.076 
12.30 
dw mull 

2.029 
4.913 
7.981 
11.23 
;iple, c( 

2.029 

4.913 

7.978 

11.10 

mnecte< 

2.029 

4.913 

7.798 

11.09 

1 domains can be 


solved. The main point to this exercise is to demonstrate how sets of coeffi- 
cients can be concantenated to form a new larger set. This new set of coefficients 
can then be solved in the same manner as before. 


Example 4.3. The following system describes the heat conduction between two 
solids. It is assumed that the heat flux and the temperatures are equal at the 
interface of the two solids. The temperature is held constant at the outer, exposed 
surfaces. Solve for the first few eigenvalues. 


d 2 u 2 
—— + A 2 u = 0 
dy 2 

d 2 v 2 

w +Xv = ° 


y=h 


= 0 
= 0 


U = V ly=0 


du 

dy 


y—0 


dv 

dy 


y—0 


(4.59) 

(4.60) 

(4.61) 

(4.62) 

(4.63) 

(4.64) 


The eigenvalues for this system can be found analytically. 



for n — 1,2, 3, . . . 


(4.65) 
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y = h 



Xj = 1 

u(y) 

y = 0 


u(Xj) 

x 1 = -1 

vfr) 

y = h 

__ 

v(Xj) 

& >< 

II 

, II 


Figure 4.1: Coordinate transformation from y to Xi and X 2 


Figure ?? shows the coordinate transformation that is needed here. 

First we will break the two layers up into two different domains. The top layer 
0 < y < h will be transformed to the first new variable — 1 < x\ < 1 and the 
bottom layer l 2 < y < 0 will be transformed into the second new variable — 1 < 
x 2 < 1. The coordinate transformation is given in the following two equations. 

*i = (£)y-i (4-66) 

= -(|)y + 1 (4.67) 

Apply the coordinate transformation to equations (4.59) through (4.64). 



d 2 u 


+ A 2 u 


dx\ 

2 d 2 v 2 

Tx | +A " 


0 

0 



(4.68) 

(4.69) 

(4.70) 

(4.71) 

(4.72) 

(4.73) 
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Two separate approximation functions are needed when applying the spectral 
method to it (xi) and v{x%). 


N 

u N (*i) = £ a n T n (x:) for - 1 < rr a < 1 (4.74) 

n=0 

M 

v m (^ 2 ) = £ b m T m (x 2 ) for — 1 < X 2 < 1 (4-75) 

771=0 

Substitute (4.74) and (4.75) into (4.68) through (4.73) and take the inner product 
of (4.68) and (4.69). 

(T) o£> + A 2 a n = 0 for n = 0, 1, . . . , AT — 2 (4.76) 

\l\ J 

0%^+A 2 fe m = 0 for m = 0, 1, . . . , M — 2 (4.77) 

£a n = 0 (4.78) 

n = 0 
Af 

£m-i)“ = 0 < 4 ' 79 ) 

m=0 

N M 

Eo.(-l) n = T.b m (4.80) 

n=0 7n=0 

(^)£« 2 an(- l) n = -(f)£mV (4.81) 

n=l ' 2 7 m =i 


If we tried to apply the matrix notation, we will run into some difficulties. The 
dependence of a n on b m and vice versa, creates difficulties in setting up the problem 
to be solved numerically. One way to work around this is to concantenate the 
set {a n : n = 0, 1, . . . , N} with the set {b m :m = 0, 1, ... , M } to create a new set, 
call it {cp : p = 0, 1, . . . ,P}. Now use the new set in (4.76) through (4.81) and 
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apply the matrix notation. 


fe) : 


E 2 


0 

I 

0 

-1) T 


(- 1 ) 


0 

2 e 2 
0 

p-N- 

-I 


V(- 1 ) B (fc) n8/ ( p-N-lfl 


Cl 

C2 


Cp-1 

Cp 


= A 2 


/ -I 

0 

0 
0 
0 

V 0 


0 \ 

-I 

0 
0 
0 
0 


Cl 

C2 


Cp-1 

cp 


J\ CP J v U U 7 V CP ) 

(4.82) 

The first row in the left and right hand side matrix of (4.82) actually has N — 
1 rows, as well as the second row. The last four rows axe exactly four rows. 
The first column is N large and the second column is M large. To clarify the 
appearance of (4.82), let N = 4 and M = 4. This gives P = 9, which maps as 
follows, {a n : n = 0, 1, . . . , 4} = {cp : p = 0, 1, . . . , 4} and {b m : m = 0, 1, . . . , 4} = 
{cp : p — 5, 7, . . . , 9}. Substitute in the values of E 2 and /, and explicitly write 
out the left hand side matrix of (4.82). 


( 0 

0 

4 ft) 2 

0 

32 GO 2 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

48 GO 2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

4(i) 2 

0 

32(f) 

0 

0 

0 

0 

0 

0 

0 

0 

24(f) 2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

48(f) 

1 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-1 

1 

-1 

1 

1 

-1 

1 

-1 

1 

-1 

-1 

-1 

-1 

-1 

1° 

Gt ) 1 


9 GO 

- 16 G0 

0 

1 

4 

9 

16 


(4,83) 
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The right hand side matrix of (4.82) looks like 

/ -1 0 0 0000 0 0 

0 -1 0 0000 0 0 

0 0 -1 0000 0 0 

0 0 0-1000 0 0 

0 0 0 000 -1 0 0 

0 0 0 0000 -1 0 

0 0 0 0000 0 -1 

000000000 
000000000 
000000000 
000000000 
\ 0 0 0 0000 0 0 


0 0 0 \ 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
-10 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 / 


(4.84) 


Again, if the two matrices were imported into a computer program, a general- 
ized eigenvalue solver could find the eigenvalues, A, as well as the solution of 
coefficients, {cp : p = 0, 1, . . . , P}. 

The following example will show how Chebyshev spectral methods can be used 
to solve for eigenvalue problems where the eigenvalue is located in the boundary 
condition. The example is a heat transfer model in a semi-infinite slab. The 
temperature is held constant at the bottom of the slab and there is Newtonian 
cooling at the top of the slab. The cooling rate is determined by the heat transfer 
coefficient, h. 


Example 4.4. Solve for the heat transfer coefficient, h, which will not give a triv- 
ial solution to the problem using Chebyshev spectral tau methods. The analytical 
solution to h is -Acot(A). 


d 2 u 

dx 2 


+ A 2 u = 0 


— — \- hu — 0 at x — 1 
dx 


u — 0 at x = -1 


(4.85) 

(4.86) 

(4.87) 


Here A is fixed and h is the eigenvalue we wish to solve for. Again we will 
substitute the Chebyshev series in for u. 

u = a n T n (x) (4.88) 

n=0 
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(4.89) 


d 2 u 
dx 2 


S °n )T n (X) 


n=0 


After substitution, the Chebyshev inner product is taken. 


N 

E (»V 


i: 


71=1 


0 for 0 < n < TV — 2 
0 

0 


(4.90) 

(4.91) 

(4.92) 


Equations (4.90) - (4.92) can now be written in matrix form. 

/ E 2 \ / 0 \ 

n 2 I I a =h —/ a (4.93) 

v (-i )•/; v ° / 

The generalized eigenvalue problem was solved using MathCad w . Table 4.2 
gives the eigenvalue h for different values of N. 

Table 4.2 The eigenvalue h for Example 4.4 for various values of N. 


N 

5 

6 

7 

8 

9 

10 

15 

Actual 

A 

8.115 

-1.885 

-0.3182 

0.4461 

0.5772 

0.5853 

0.5883 

0.5883 


Note that in this example, there was only one eigenvalue to be found. In 
general, the number of eigenvalues that will be found will be equal to the rank 
of the left hand matrix or the right hand matrix, whichever is smaller. Here the 
rank of the right hand matrix is equal to 1, which is the number of eigenvalues 
found. 

The next and last example will demonstrate an actual application of this 
method to the Nield problem [7]. The Nield problem is the solution of the onset 
of fluid convection in a single layer of liquid. The liquid layer is either heated 
from above or below by a heating plate. The upper surface of the liquid is free to 
deflect and is exposed to air (or some other gas). 

Example 4.5. Given the following system, solve for the critical Marangoni num- 
ber (Ma) which corresponds with the onset of fluid convection. 

(D 2 - u 2 f W = u 2 Rae 
(D 2 - uj 2 ) 9 = -W 
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(4.94) 

(4.95) 



The system needs six boundary conditions plus an additional equation for the 


deflection of the gas-liquid interface, C- 

W(- 1) = 0 (4.96) 

DW(-l) = 0 (4.97) 

©(-1) = 0 (4.98) 

W(0) = 0 (4.99) 

(D 2 - 3u> 2 ) DW (0) - Ra + - ^ w 2 ? = 0 (4.100) 

D 2 W (0) - u?Ma [q (0) - 0 (0)] = 0 (4.101) 

D© (0) + L [0 (0) - <r] = 0 (4.102) 


The system of equations has been non-dimensionalized and linearized using a Un- 
ear stability analysis. The time derivative has been ehminated using a Fourier 
expansion. Because we are only interested in the steady solutions, the time con- 
stant is set equal to zero. The problem will be set up so that the Marangoni 

d _ 

Number, Ma, is the eigenvalue. Here D is the derivative — , W is the velocity, 0 

is the temperature, u) is the wave number corresponding to the perturbation, Ra 
is the Rayleigh number, G is the Weber number, C is the Crispation number, Ma 
is the Marangoni number, L is the Biot number and c is the surface deflection 
term. For a complete derivation, the interested reader is referred to the paper by 
Nield. In Nield’s paper, the surface deflection is neglected. Here the surface is 
allowed to deflect. Values for the deflecting case have been calculated [10]. 

Again we need to stretch the coordinates from 0<z<lto— 1 < x < 1, 
using the stretching function, x = 2z + 1. Apply the Chebyshev polynomial 
interpolation. 

W = f> n r n (x) (4.103) 

n=0 

M 

© = E b m T m (x) (4.104) 

m=0 

After substitution, the inner product is taken. 

(l6E 4 — 8o; 2 E 2 + u; 4 l) ■ a. — u 2 Ral ■ b = 0 (4.105) 

(4E 2 - u; 2 l) b + a • I = 0 (4.106) 
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N 


EMK = 0 

(4.107) 

n=0 

E(-l )*”«’«. = 0 

(4.108) 

71=1 

M 

= 0 

(4.109) 

o 

II 

cf 

*W1 

O 

II 

6 

(4.110) 


N [ n 2 


8E 

n=3 L 


- (n 2 - 1) (n 2 - 4) 


Rcl + 


N 


— 6u 2 ^2 n2(1 n 
n= 1 

G + UI 2 1 2 

= 0 


JV „2 / M \ 

4^2 — (n 2 - l) a n - u 2 Ma ( jn <r - ) 

n= 2 J \m=0 / 

Af M 

2 m 2 6 m + L £ 6 m - Lq 


0 

0 


(4.111) 

(4.112) 

(4.113) 


m=0 


E 4 can be found by using equation (3.49). The strange products in (4.111) and 
(4.112) are the result of the third and second order derivative in equation (3.35). 

There are several interesting points to be made at this time. In (4.105), the in- 
ner product of the equation restricts the number of terms that b = {6 m : m = 0, 1, . 
could have, which is the number of inner products taken. Here we will choose 
N — 4. The number N — 4 corresponds to the choice of N and that there are 
4 boundary conditions associated with (4.105). The same situation holds for 
(4.106). Again we will choose N but here there are only two boundary conditions 
associated with the equation, which gives N — 2 equations. Equation (4.105) ap- 
plies restrictions to the choice of M. When the inner product of the equation is 
taken, @/v must lie in Wjv and vice versa. 

For simplification and in general, we will let TV = M and truncate the in- 
dices accordingly. With this in mind, equations (4.105) through (4.112) can be 
rewritten. 


(l6E 4 -8a; 2 E 2 +u; 4 l) a-u 2 Ral b = 0 (4.114) 

(4E 2 -w 2 l)b + a-I = 0 (4.115) 
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0 


(4.116) 


N 


rl— 0 

\n+l „2 


^(-ir^V = 0 


71= 1 


JV 


£(-1)6™ = 0 
m=0 

AT 

5>n = 0 


71=0 


N 

8E 

71=3 


" n 2 


nr 

27 


(n 4 — 5 n 2 + 4) 


N 


Ra + 


— 6 U> 2 ^2 n2a r> 

71= 1 

G + u 2 ' 


uj 2 q = 0 


N 


4 — (n 2 - l) a n - uj 2 Ma 51 * “ h ™ ) = 0 

n= 2 6 \m=0 / 

Now write (4.113) through (4.120) in matrix form. 

—uj 2 Ra I 
4E 2 - u 2 l 
0 
0 


/ 16E 4 - 8u; 2 E 2 + u; 4 I 

I 

(-1 )"/ 

( — l) n n 2 1 
0 

I 

8S| (n 4 — 5n 2 + 4) — 6 u 2 n 2 


V 


(n 2 - 1 ) 

0 


(- 1 )” 

0 

0 

0 

(2n 2 + L) I 


0 

0 

0 

0 

0 

0 


\ 


-(Ra + e^-)u 2 
0 

-L 


) 


/ a 0 
a i 

On 

bo 




Ma 


( 0 0 

0 0 

0 0 
0 0 
0 0 
0 0 
0 0 


0 \ 
0 
0 
0 
0 
0 
0 


Oo 
a i 


&71 

bo 


0 —uj 2 I u 2 

Vo o o ) 


bn 

\ * 


(4.117) 

(4.118) 

(4.119) 

(4.120) 

(4.121) 

\ 

(4.122) 
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The first row of both matrices in (4.122) has actually N — 3 entries and the second 
row has N — 1 entries. The first and second column has N entries, and the last 
column has only one entry. The FORTRAN program used to set up and solve for 
the above generalized eigenvalue problem is given in the appendix. One of the 
subroutines was called from an IMSL subroutine library. 

5. Summary 

I hope that this paper has helped to demonstrate the usefulness of Chebyshev 
Spectral Tau methods. Of course there is no such thing as a free lunch. It 
becomes obvious that the implementation of this method is more complicated than 
a finite difference method. The Chebyshev Spectral Tau method is most useful 
for problems which may need to be solved for many times or when a parameter in 
the problem needs to be changed often. The usefulness of spectral methods was 
also demonstrated for eigenvalue problems. 

In section two, the Chebyshev polynomial properties and the methods of sum- 
mation were given as background material for the third section. Section three de- 
veloped the recurrence relationships between functions expanded into Chebyshev 
polynomials and the derivatives of those expansions. The relationships between 
the derivatives of the expansion coefficients were used in section four in the spec- 
tral tau technique. These recurrence relationships can be found in the appendix 
of Gottlieb and Orszag’s book [5] for various linear operators. If the relationship 
you are interested in is not tabulated, the techniques given in this paper can be 
used to generate them. The appendix gives a copy of the FORTRAN program 
used to solve for the eigenvalues of Nield’s problem in Example 4.5. Some of the 

newer higher level software could also be used to solve these equations. Good 

(r) (r) (r) 

examples of these programs are MathCad , Maple w and Matlab . 

The intent of this report is to help the beginner or uninitiated start using 
Chebyshev spectral methods. There are several aspects which were not dealt 
with. Some of these are singularities in the governing equations, efficient solver 
methods, such as fast Fourier transforms, and stability. These topics are covered 
in Spectral Methods in Fluid Dynamics [2]. It appears that fourth order and higher 
differential equations will become unstable. This problem can be worked around 
by splitting fourth order equations into two second order equations. Even though 
the fourth order equations are unstable, most problems can still be solved using 
the Chebyshev Spectral Tau technique. A good example is given in Numerical 
Analysis of Spectral Methods [5] on page 144-145. Some higher order problems, 
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however do exhibit unstable behavior. I have discovered a system which gives 
incorrect eigenvalues for only certain values of the system parameters. This system 
was similar to the Nield problem but dealt with two liquid layers with deflecting 
surfaces. Whenever possible, it is advised to deal with second order equations 
only. 
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6. Appendix 

PROGRAM NIELD 

INTEGER I, J, N, M, BC, II 

REAL PI, OMEGA, OMEG2, OMEG4 

PARAMETER (N = 20, M = 2*N+1, PI = 3.141592654) 

REAL B2(N, N), B4(N, N), A(M,M), B(M,M), BETA(N) 

REAL VECON(N), VECl(N), VECIN(N), VEC2(N), VEC3N(N) 
REAL L, MA, RA, OMEGA 
REAL TNI, TN2, TN3, TN4, TNO 
COMPLEX ALPHA(M), EVEC(M) 

* Define the matrices used in defining the large matrices ********** 
*** INTMAT initializes the matrix entries to zero. *************** 
CALL INTMAT(B2, 1, N, 1, N) 

CALL INTMAT(B4, 1, N, 1, N) 

CALL INTMAT(A, 1, M, 1, M) 

CALL INTMAT(B, 1, M, 1, M) 
jjc afc i^ "|^ g iv e s the matrix 
CALL CHDER2(B2, N, N) 

CALL CHDER4(B4, N, N) 

***CHDERB gives the boundary condition vector ************** 
CALL CHDERB ( VECON , 0, N, -1) 

CALL CHDERB (VEC1, 1, N, 1) 

CALL CHDERB(VEC1N, 1, N, -1) 

CALL CHDERB(VEC2, 2, N, 1) 

CALL CHDERB(VEC3, 3, N, 1) 

^Define the parameters to be used. OMEG2, OMEG4, TNI through * 
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*TNP1 are used for optimization puposes *********************** 
OMEGA = 2.13 
L = 1.04 
G = 0.452 
C = 50.9E-5 

OMEG2 = OMEGA * OMEGA 

OMEG4 = OMEG2 * OMEG2 

TNI = 2 * N - 1 

TN2 = 2 * N - 2 

TN3 = 2 * N - 3 

TN4 = 2 * N - 4 

TNO = 2 * N 

** Row 1 Column 1 ******* 

DO 20 I = 1, N-4 
DO 10 J = 1+1, N 

A(I,J) = 16.0 * B4(I,J) - 8.0 * OMEG2 * B2(I,J) 

10 CONTINUE 
A(I,I) = OMEG4 

20 CONTINUE 

** Row 1 Column 2 ******* 

DO 30 I = 1, N-4 
A(I, I+N) = -OMEG2 * RA 
30 CONTINUE 
** Row 2 Column 1 ******* 

DO 40 I = 1, N-2 
A(I+N-4,I) = 1.0 
40 CONTINUE 
** Row 2 Column 2 ******* 

DO 60 I = 1, N-2 

11 = I+N-4 

DO 50 J = 1+1, N 
A(II,J+N) = 4.0 * B2(I,J) 

50 CONTINUE 
A(II, I+N) = -OMEG2 
60 CONTINUE 

** Boundary Conditions ********** 

DO 100 J = 1, N 
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A(TN4,J) = VECON(J) 

A(TN3,J) = VECIN(J) 

A(TN2,J+N) = VECON(J) 

A(TN1,J) = 8.0 * VEC3(J) - 6.0 * 0MEG2 * VECl(J) 

A(TNO,J) = 4.0 * VEC2(J) 

A(M,J+N) = L 
A(M,J) = 2.0 * VECl(J) 

100 CONTINUE 

A(TN1,M) = -OMEG2 * (RA + (G + OMEG2) / C) 

A(M,M) = -1.0 

********* * ***** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

*** Now assign values to the B matrix ***** 

DO 200 J = 1, N 
B(TN0,J+N) = -OMEG2 
200 CONTINUE 
B(TN0,M) = OMEG2 

***************************************************************** 

*The subroutine GVLRG is linked and called from ISML library. 

*M is the size of the matrix, ALPHA is a vector of numerators of the 

^eigenvalues, and BETA is a vector of denominators of the eigenvalues. 
***************************************************************** 

CALL GVLRG (M, A, M, B, M, ALPHA, BETA) 

DO 200 J = 1, M 
IF (BETA(J) .EQ. 0) THEN 
EVEC(J) = 12345 
ELSE 

EVEC(J) = ALPHA(J) / BETA(J) 

END IF 

200 CONTINUE 
END 

*******************************************************„.*.,,.****,,;,,<* 

* This subroutine evaluates the derivatives of a Chebyshev polynomial 

* at a boundary point of 1 or -1. It returns the vector of 

* coefficients that are associated with it in VECTOR. 

* VECTOR is the vector returned 

* ORDER is the order of the derivative 

* LENGTH is the length of the vector VECTOR 


43 



* EVAL is the point at which the derivative is evaluated 
SUBROUTINE CHDERB (VECTOR, ORDER, LENGTH, EVAL) 
INTEGER ORDER, LENGTH, K, N, EVAL, N2 
REAL VECTOR(LENGTH), PROD 
IF (EVAL .EQ. -1) THEN 
DO 20 N = 0, LENGTH-1 
PROD = 1.0 
N2 = N * N 

DO 10 K = 0, ORDER - 1 

PROD = PROD * (N2 - K*K) / (2.0 * K + 1) 

10 CONTINUE 

VECTOR (N+l) = (-1.0)**(N+ORDER) * PROD 
20 CONTINUE 

ELSE IF (EVAL .EQ. 1) THEN 
DO 40 N = 0, LENGTH-1 
PROD = 1.0 
N2 = N * N 

DO 30 K = 0, ORDER - 1 
PROD = PROD * (N2 - K*K) / (2.0 * K + 1) 

30 CONTINUE 
VECTOR(N-hl) = PROD 
40 CONTINUE 
ELSE 

PRINT *, ’Boundary condition evaluated at’ 

PRINT *, ’value other than 1 or -1.’ 

PAUSE ’Executing in subroutine CHDERB’ 

END IF 
END 

***************************************************************** 

* This subroutine finds the coefficient matrix associated with 

* the first derivative of the Chebyshev series 

* ROWS is the number of rows in the matrix 

* COLUMNS is the number of columns in the matrix 

* I and J are indices used. 

* MATRIX is the matrix that is passed back to the main program 

* MATRIX has size (l:ROWS, T.COLUMNS) 
***************************************************************** 
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SUBROUTINE CHDER1 (MATRIX, ROWS, COLUMNS) 

INTEGER I, K, ROWS, COLUMNS 
REAL MATRIX ( ROW S , COLUMNS) 

CK = 2 

DO 10 I = 1 , ROWS 
DO 20 J = I + 1 , COLUMNS , 2 
MATRIX (I, J) = 2 / CK * (J - 1) 

20 CONTINUE 
CK = 1 

10 CONTINUE 
END 

********** ****************** ********* ******* ********************* 

* This subroutine finds the coefficient matrix associated with 

* the second order derivatives of the Chebyshev series 

* ROWS is the number of rows in the matrix 

* COLUMNS is the number of columns in the matrix 

* I and J are indices used. 

* MATRIX is the matrix that is passed back to the main program 

* MATRIX has size (l:ROWS, lrCOLUMNS) 

* * * * * * * * * * * * 5k * * * * * * * * * * * * * * * * 5k * * 5k * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

SUBROUTINE CHDER2 (MATRIX, ROWS, COLUMNS) 

INTEGER I, K, ROWS, COLUMNS 
REAL MATRIX(ROWS, COLUMNS) 

CK = 2 

DO 10 I = 1 , ROWS 
DO 20 J = I + 2 , COLUMNS , 2 
MATRIX (I, J) = ((J-l) * ((J-l)**2 - (I-l)**2)) / CK 
20 CONTINUE 
CK = 1 

10 CONTINUE 
END 

***************************************************************** 

* This subroutine finds the coefficient matrix associated with 

* the fourth order derivatives of the Chebyshev series 

* MATRIX has size (IrROWS, l:COLUMNS) 

* MATRIX = the matrix to be passed back assigned with the value 

* of fourth order coefficients 
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* ROWS = the length of the rows of MATRIX 

* COLUMNS = the number of the columns of MATRIX 
***************************************************************** 

SUBROUTINE CHDER4(MATRIX, ROWS, COLUMNS) 

INTEGER I, J, ROWS, COLUMNS 

REAL MATRIX(ROWS, COLUMNS), A, B, C, D 

CK = 2 

DO 20 I = 0, ROWS-1 
D = 1**2 * (1**2 - 4)**2 
DO 10 J = I + 4, COLUMNS-1, 2 
A = J**2 * (J**2 - 4)**2 
B = 3.0 * 1**2 * J**4 
C = 3.0 * 1**4 * J**2 

MATRIX(I+1, J+l) = J*(A - B + C - D) / (CK * 24) 

10 CONTINUE 
CK = 1 

20 CONTINUE 
END 

***************************************************************** 

* This subroutine initializes all the entries of a matrix 

* ROWS is the number of rows in the matrix 

* COLUMNS is the number of columns in the matrix 

* I and J axe indices used. 

***************************************************************** 
SUBROUTINE INTMAT (MATRIX, ROW1, ROWEND, COL1, COLEND) 
INTEGER I, J, ROW1, ROWEND, COL1, COLEND 
REAL MATRIX ( ROW 1 : ROWEND , COLl:COLEND) 

DO 20 J = COL1 , COLEND 
DO 10 I = ROW1 , ROWEND 
MATRIX (I, J) = 0.0 
10 CONTINUE 
20 CONTINUE 
END 
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